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Abstract 

We study the collective behaviour of an ensemble of coupled motile elements whose 
interactions depend on time and are alternatively attractive or repulsive. The evo- 
lution of interactions is driven by individual internal variables with autonomous 
dynamics. The system exhibits different dynamical regimes, with various forms of 
collective organization, controlled by the range of interactions and the dispersion 
of time scales in the evolution of the internal variables. In the limit of large inter- 
action ranges, it reduces to an ensemble of coupled identical phase oscillators and, 
to some extent, admits to be treated analytically. We find and characterize a tran- 
sition between ordered and disordered states, mediated by a regime of dynamical 
clustering. 
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1 Introduction 



Large ensembles of coupled dynamical systems are a basic tool for the study 
of emergent collective behaviour in systems of interacting agents [1]. Self- 
organization in spatially distributed ensembles is revealed by the appearance 
of spatiotemporal structures, which may exhibit complex evolution. Instead, 
coupled ensembles for which space is not relevant to the dynamical laws be- 
come organized in temporal patterns, typically, in some form of synchronous 
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evolution. Much attention has been focused on synchronization phenomena 
over the last years [2]. They have been identified in a broad class of natural 
and artificial systems, and many of them have been successfully explained by 
means of relatively simple models. 

The degree of collective organization in a coupled ensemble, measured by the 
correlation between individual motions, varies considerably between different 
kinds of synchronization patterns. Most studies on these phenomena have 
dealt with patterns of highly correlated evolution, such as full and phase syn- 
chronization. Indeed, these forms of synchronization are essential to the func- 
tioning of some artificial systems, and have been observed in certain insect 
populations [3,4]. On the other hand, they are expected to play a less relevant 
role in most biological systems, where the complexity of collective functions 
requires a delicate balance between coherence and diversity [1]. Consider, for 
instance, the brain, where highly coherent activity patterns are only realized 
under pathological states, such as during epileptic seizures. Many biological 
systems consisting of interacting agents, ranging from biomolecular complexes 
to social populations, are normally found in configurations where the ensemble 
is segregated into groups with specific functions. While the evolution of indi- 
vidual elements is highly correlated inside each group, the collective dynamics 
of different groups is much more independent. Clustering has been modeled by 
means of coupled ensembles with specific individual dynamics or interactions 
[5,6,7,8,9]. 

Usually, clustering is a dynamical process, where groups may preserve their 
identity in spite of the fact that single elements are continuously migrating 
between them. Intermittent formation of clusters and migration of elements 
has been observed under the action of noise [10,11,12]. The individual motion 
towards or away from clusters may also be controlled by the internal state 
of each element, which favors or inhibits grouping with other elements. This 
is observed in natural phenomena ranging from complex chemical reactions, 
where biomolecules react with each other only when they have reached appro- 
priate internal configurations [13], to social systems, where the appearance of 
organizational structures requires compatibility between the individual chang- 
ing states of the involved agents [14]. Focusing on this motivation, in this paper 
we study a model of coupled motile elements where interactions depend on 
the internal state of each element. Individual internal variables change with 
time according to a prescribed autonomous evolution law, so that interactions 
are also time-dependent, and result to be alternatively attractive or repulsive. 
This intricate interaction pattern induces several complex regimes of collective 
evolution, including dynamical clustering. Clustering in systems of interact- 
ing moving particles with internal dynamics has been demonstrated when the 
internal states of different elements are mutually coupled and their individual 
evolution is chaotic [15] . In such case, the synchronization of internal variables 
results into spatial organization, through the effect of those variables on the 
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degrees of freedom associated with motion. Here, we show that similar forms of 
collective behaviour are already possible with simple, autonomous oscillatory 
dynamics for the internal variables. 

In the next section, we introduce the model and present numerical results 
for the case of two-dimensional motion in a bounded domain. The occurrence 
of different dynamical regimes is controlled by the range of interactions and 
by the dispersion in the time scales associated with the internal dynamics. 
In Section 3, the limit of large interaction ranges is analyzed. In this limit, 
the system reduces to an ensemble of identical phase oscillators with internal 
states. As the dispersion in the time scales of internal variables is varied, it 
exhibits an order-disorder transition between a configuration where the phase 
and the internal variable are highly correlated and a uncorrelated state. This 
critical transition, which we are able to describe analytically, is mediated by a 
regime of dynamical clustering. Both the transition and the clustering regime 
are characterized by means of suitably defined order parameters. In the final 
section, we discuss our results with emphasis in the features that might be 
also observed in other systems of similar kinds. 



2 Dynamical clustering of motile elements 

2.1 Model 

We consider an ensemble of N motile elements at positions rj (i = 1, . . . , A^). 
The elements interact through pair forces F(rj— r^) depending on their relative 
positions. Moreover, the interaction of each pair is modulated by a function 
V{0i,0j) of the variables 9i and dj, which characterize the internal state of 
the respective elements. This modulation aims at describing the effects of the 
internal states on the spatial dynamics. The motion of the elements is assumed 
to be overdamped, so that the equation of motion for element i is 



From the viewpoint of the numerical resolution of Eqs. (1), a convenient choice 
for the interactions a truncated elastic force. 



r,; = 





k{ri - Tj) if \ri - rj\ < ro. 



(2) 







otherwise. 



Here, k > is the elastic constant and Tq is the interaction range. 



3 



As for the internal state described by the variable 6i, we assume that each 
element behaves as an autonomous oscillator with constant frequency fli, much 
like a "biological clock." Biological periodic phenomena at several levels, from 
intracellular chemical reactions to complex long-period rhythms, are in fact 
well described by such elementary cychc processes [16]. The internal variable 
represents thus a phase, and evolves in time as 

ei(t)^nit + ei(o). (3) 



To take into account inhomogeneity in the ensemble, the frequencies fli are 
drawn from a given distribution g{fl). Here, we take the Gaussian 



v27rcr^ 



2(72 



(4) 



The autonomous dynamics of the internal phases 9i represents cyclic processes 
which, in our model, are not influenced by external factors but affect the 
interaction between elements. Specifically, the interaction weight V{9i, 9j) in 
Eq. (1) is chosen as 

viOi, Oj) = cosiOi - Oj) = cos[{ni - nj)t + ^°.], (5) 



with — 9i{0) — 9j{0). With this choice, the interaction between elements i 
and j changes its sign periodically, in a time scale of order (Qj — Q^)"^. Thus, it 
is alternatively attractive or repulsive, respectively favoring or inhibiting the 
aggregation of elements. At the level of the whole ensemble, this variations 
give rise to a complex time-dependent interaction pattern which, as shown 
below, induces a variety of dynamical regimes, including dynamical clustering. 
Note that, according to Eq. (5), interactions between pairs of elements do 
not depend separately on the individual frequencies, but on the frequency 
differences — flj. Thus, in Eq. (4), we can fix Qq ~ 0. 

As a final ingredient of the model, it is necessary to consider that the ensemble 
is confined within a bounded spatial domain. In fact, when mutually attracting 
elements have aggregated into clusters, we expect that the remnant repulsion 
between elements in different neighbor clusters drives them away from each 
other. In the long run, this effect would lead to the total dispersion of the 
ensemble. Therefore, we assume that the system evolves into a closed domain 
with refiecting boundaries. 
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2.2 Numerical results 



We have extensively studied the above model, by solving Eqs. (1) numerically 
for a two-dimensional system, where an ensemble oi N — 100 elements is 
confined within a circle of radius rmax- Simultaneous rescaling of frequencies 
and time makes it possible to fix the elastic constant to /c = 1. Moreover, we 
choose Tmax as the unit of distance. In this way, the only relevant parameters 
in the model are the interaction range Tq and the frequency dispersion a. 

Truncation errors in the numerical resolution of the equations of motion (1) 
may cause that, under the effect of attractive interactions, two elements col- 
lapse to numerically indistinguishable positions. This is expected to happen, 
in particular, when the evolution of internal phases is slow (small a) so that 
attractive interactions can act during sufficiently long times. Under such con- 
ditions, in view of Eq. (2), the mutual interaction of the collapsed elements will 
vanish, and the forces applied on each of them by any other element will be 
identical. Consequently, the collapsed elements will remain "stuck together" 
for the remaining of the evolution. To avoid this numerical artifact, we intro- 
duce in Eqs. (1) a small additive-noise term which prevents the spatial collapse 
of elements and, otherwise, does not affect the dynamics. 

Figure 1 shows four successive snapshots of the ensemble, in a single realization 
with To = 0.3 and a = 0.1. We find that many elements are entrained in several 
compact clusters, while others are distributed in space along curved lines, 
forming "strings." These "strings" result from the instabilization of clusters, 
as illustrated in Figs, la and b. Generally, the distribution of internal phases 
6i of the elements belonging a given chister spans a total angle less than tt, but 
this total angle is often larger than 7r/2. This implies that, inside a cluster, 
the difference of internal phases between a pair of elements can be larger than 
7r/2, and their interaction can thus be repulsive. However, the two elements 
can be part of the same cluster due to the presence of other elements with 
intermediate internal phases. These elements mediate the interaction between 
the otherwise repulsive pairs, in such a way that the effective force is attractive 
towards the cluster. 

Once a cluster is formed, and after a certain transient has elapsed, the total 
angle spanned by the distribution of internal phases 6i grows with time, due to 
the individual evolution of each phase, Eq. (3). This growth takes place within 
time scales of order min^ |Qj — Under some simplifying assumptions, it 

can be shown that a cluster becomes unstable when the total angle spanned by 
the internal phases reaches vr. Suppose that, at a given moment, the internal 
phases of the n elements belonging to a cluster are uniformly distributed over 
the interval (0,0). Suppose also that all the elements are at the same point 
in space, except for a test element at a small distance Sr (with \6r\ < tq) from 
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Fig. 1. Four snapshots of an ensemble of size N = 100, with tq = 0.3 and a = 0.1, at 
times (a) t = 5100, (b) t = 5110, (c) t = 5120, and (d) t = 5130. Each dot represents 
the position rj = (xj, yi) of an element on the plane, while the orientation of the 
stick starting at each dot indicates the value of the corresponding internal phase 6i. 
Circles show the system boundary. All the elements of the cluster marked with an 
arrow in (a) are part of the "string" marked in (b). 

the cluster. If n is sufficiently large, taking into account Eqs. (1) to (5), the 
equation of motion of the test element can be written as 



5r ——kSr'Y^ cos{9i — 9o) ^ knSr — I cos(^o — 0)d9 — 
^ " 

^^^ sin(go-e)-singo ^^^ (6) 

where 6q is the internal phase of the test element. The factor multiplying Sr in 
the right-hand side of Eq. (6) is negative for all Oq G (0, 0) if < < vr. Under 
these conditions, the force on the test element is attractive. As reaches vr, 
however, the factor of 6r becomes positive at 6*0 = and n. This implies that, 
as soon as > tt, the elements with extreme internal phases, 9 ^ or 9 Ki n, 
will be repelled from the cluster. Due to their mutual repulsion, elements with 
internal phases close to ^ = and 9 = n will move away from the cluster in 
opposite directions. Moreover, as the distribution of phases in the remaining of 
the cluster keeps broadening and new elements reach the stability limits, they 
will be attracted by those elements with similar phases which have just left 
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Fig. 2. Density plot of the histogram of pair distances as a function of time, for the 
same reahzation as in Fig. 1 {N = 100, tq = 0.3, a = 0.1), along 200 time units. 
Darker tones correspond to higher concentrations. 

the cluster, and therefore move in their direction. This mechanism gives rise to 
the "strings" of elements when a cluster becomes unstable and disintegrates. 

The dynamics of clustering is effectively illustrated by the distribution of dis- 
tances between elements [6,7]. Calculating the N[N — l)/2 Euclidean pair 
distances 



dij 



(7) 



for all i and j > i, we construct a histogram where the height of the column 
with base {d, d + Ad) is proportional to the number of pair distances within 
that interval. Compact clusters cause the appearance of sharp peaks in the 
histogram, at ^ and other distances, while non-clustered elements give rise 
to a smooth background. Figure 2 shows, as a density plot, the evolution of 
a 200-column histogram of pair distances d e (0, 2) {Ad — 0.01) for the same 
numerical reahzation as in Fig. 1, between t — 9800 and 10000. As expected, 
the region d < ro is strongly depleted (except for the peak at ^ 0), because 
pairs of elements at such distances either become entrained into clusters or are 
mutually repelling. For d > Tq, dark horizontal lines reveal localized compact 
clusters. As time elapses, these clusters can become less well-defined, and they 
are seen to move, disperse, coalesce, or split to form smaller groups. 

This complex dynamics, driven by the autonomous evolution of the internal 
phases, calls for a statistical description in terms of a characterization of the 
partition of the ensemble into clusters. With this aim, we introduce two order 
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parameters [6]. The parameter R is the time-avcraged fraction of pair distances 
smaller than a certain small threshold 8. It reads 



where (■) indicates average over long times, and H{z) is the step function: 
H{z) — 1 for z > 0, and H{z) — otherwise. The second parameter, S, is 
given by the timc-averaged fraction of elements which have at least another 
element as a distance smaller than 5. It can be calculated as 



Associating the threshold 5 with the (maximal) size of a cluster, both param- 
eters are statistical measures of the degree of clustering in the ensemble. Their 
values are close to zero if most elements are not entrained, and reach unity 
if all elements belong to clusters. They are, however, independent quantities. 
For a fixed value of 5, which measures the total population inside clusters, 
R is small if the number of clusters is large, and vice versa. In other words, 
R characterizes the average number of elements per cluster. An approximate 
evaluation of the number of clusters M is obtained by assuming that all clus- 
ters have the same population, as M = S"^ / R [6]. 

To characterize clustering in our system, we have numerically calculated the 
order parameters R and S as functions of the frequency dispersion a and the 
interaction range Tq in a system of iV = 100 elements, taking 5 = 0.01. Figure 
3 shows R and S versus a for tq = 0.3. For this intermediate value of the 
interaction range, the two order parameters display a moderate decrease as 
the frequency dispersion grows. As frequency differences become larger, the 
typical length of the time intervals during which the interaction of a given pair 
of elements is attractive shortens. Consequently, shorter times are available 
for the aggregation of elements into clusters and, thus, the average degree of 
clustering diminishes. The approximate number of clusters, calculated as M = 
S'^/R, is practically constant, M 20, in the considered interval of variation 
for a. The same behaviour is found for smaller values of the interaction range. 
On the other hand, for large values of Tq, the dependence with the frequency 
dispersion becomes less trivial. We analyze this specific situation in the next 
section. 

The dependence of the clustering order parameters R and S on the interaction 
range Tq for fixed a is more interesting, as shown in Fig. 4 for a = 0.1. We find 
that both parameters display a maximum, at Tq ~ 0.27 for R and Tq ~ 0.23 for 
S. For smaller values of the interaction range, the two parameters rapidly drop 




(8) 




(9) 
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Fig. 3. The order parameters R and S as functions of the frequency dispersion a, 
for an interaction range ro = 0.3. 
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Fig. 4. The order parameters R and S as functions of the interaction range ro, for 
a frequency dispersion cr = 0.1. 

to zero, while for larger tq they exhibit a smoother decay. The approximate 
number of clusters, M = S'^/R, grows steadily from M 1 for tq ~ and 
attains a sharp maximum, M ^ 30, at ro ~ 0.25. For larger interaction ranges, 
it soon reaches an almost constant value, M 20. 



In order to explain the non-monotonic dependence of the clustering order 
parameters with the interaction range, it is useful to inspect the instantaneous 
spatial distribution of the ensemble for different values of ro (Fig. 5). For 
sufficiently small ro, there is little chance that two elements are simultaneously 
found within the interaction range. In our two-dimensional system, the average 
distance d between elements can be evaluated as d — A/N, where A is the 
area of the domain accessible to the ensemble. Taking A — irr^^^ and rmax = 1, 
we find d ~ 0.18. Figures 5a and b show that for < d only a few elements 
form clusters, whereas for Tq > d most of the ensemble is clustered in many 
groups. For such values of the interaction range, consequently, the clustering 
order parameters increase with ro (Fig. 4). 

As the interaction range grows, the collective organization of elements in larger 
spatial structures becomes possible. In particular, the "strings" associated 
with unstable clusters, discussed above, entrain now a substantial part of the 
ensemble (Fig. 5c). The contribution of these structures to the clustering or- 
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Fig. 5. Four snapshots of an ensemble of size N = 100 at times t = 10000, with 
a = 0.1 and different interaction ranges: (a) vq = 0.1, (b) ro = 0.2, (c) ro = 0.3, 
and (d) ro = 0.5. 

der parameters is expected to be considerably weaker than that of compact 
clusters, which explains the dechne of R and S for ro ^ 0.3. For even larger 
interaction ranges the dynamics is dominated by the remnant repulsive inter- 
action between compact clusters. When Tq is sufficiently large, all the elements 
are found most of the time at the boundary of the spatial domain (Fig. 5d). 
Note that elements are still organized in space according to their internal state: 
generally, neighbouring elements have similar internal phases. This form of or- 
ganization in the limit of large interaction ranges is treated in detail in the 
next section. 

The change in the collective spatial distribution as the interaction range grows 
is also apparent in the distribution of pair distances dij. Figure 6 shows 200- 
column normalized histograms of pair distances for the snapshots shown in 
Fig. 5. When the degree of clustering is small (Fig. 6a) the distribution of 
pair distances is rather smooth, and no significant contributions are present at 
d ~ 0. As the interaction range grows and clustering becomes more important, 
the distribution spreads and develops sharp spikes. Note, in particular, the 
high peak at d ^ (Fig. 6b and c). For even larger Tq the degree of clustering 
recedes, the peak at c? ~ becomes smaller, and other spikes are replaced 
by broader structures (Fig. 6d). The highest peak, centered at o? ~ 1.8, is 
a direct byproduct of the concentration of elements in the boundary of the 
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Fig. 6. Normalized histograms of pair distances for the snapshots of Fig. 5. (a) 
ro = 0.1, (b) ro = 0.2, (c) ro = 0.3, and (d) ro = 0.5. 

spatial domain. 

The fact that, in the hmit of large interaction range, the ensemble is mainly 
concentrated in the boundary of the accessible domain suggest a reduced rep- 
resentation of the problem, where the position of each element is given by a 
linear coordinate along the boundary. In the case of the circular domain con- 
sidered above, such approximation yields a set of dynamical equations which 
are formally equivalent to the model equations of coupled phase oscillators. As 
shown in the next section, this reduced representation admits to be studied 
analytically to a considerable extent. 



3 Phase oscillators with time-dependent coupling 

3.1 Reduction to phase variables 

When the range of interactions is larger than the linear size of the system, all 
the elements are subject to the action of the whole ensemble at all times. As 
discussed in the preceding section, the regime of large ro is characterized by the 
repelling interaction of clusters, which drives the elements to the boundaries of 
the accessible spatial domain. For the case of the circular domain considered 
above, the position of an element at the boundary can be characterized by the 
angle 0j defined by the element and an arbitrarily chosen origin, with vertex 
at the center of the circle. Assuming that the elements do not abandon the 



11 



Fig. 7. Reduction of the circular-domain problem to phase variables, in the limit of 
large interaction range. 

boundary at any time, the evolution of 4>i is governed by 

<i)i = -Xi sin (f)i + yi cos 0i = Vi^i, 0j)Fij, (10) 



where F-j is the projection of the pair force F(rj — r^) parallel the circular 
boundary [cf. Eq. (1)]. Figure 7 shows the geometry of the problem. With 
ro > 2rmax = 2, the modulus of the force (2) between any pair of elements is 
given by 

Fij = k dij = A;^2[l - cos(0,- - 00] = sin (^^^^^ ■ (H) 



The angle between the force and the tangent to the circle at the position of 
element i is {(pj — 0j)/2. Therefore, 

= Fij cos ( ) = A;sin(0,- - 0,)- (12) 



The equation of motion for 0j becomes 

^i = kY,v{ei,ej)sm{<pj-<p,). (13) 



Remarkably, this equation is formally equivalent to the model equation for 
an ensemble of phase oscillators with identical natural frequencies a; = 0, 
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where the state of each oscihator is characterized by its phase ^j. The col- 
lective behaviour of ensembles of phase oscillators has been thoroughly an- 
alyzed during the last two decades. With V{9i,9j) = 1 for all i,j, Eq. (13) 
reduces to the Kuramoto model for globally coupled identical oscillators [17]. 
In this case, full synchronization {(f)i{t) — (f)j{t) for all i,j) is a stable asymp- 
totic state for any positive value of k. Time-independent interaction weights, 
where V{6i,9j) = Vij is constant but can have different values for each oscil- 
lator pair, have been shown to induce typical features of disordered systems, 
such as glassy-like behaviour, frustration, and algebraic relaxation towards 
equilibrium [18,19,20]. As we show in the following, when the interaction be- 
tween phase oscillators is modulated by the time-dependent interaction weight 
V{6i,6j) defined by Eqs. (3) and (5), various forms of collective evolution - 
including dynamical clustering- emerge. The transition between these regimes 
is controlled by a, the frequency dispersion of the internal phases di, which is 
now the only parameter left in the model. 

3.2 Order- disorder critical transition 

We begin our analysis of Eq. (13) by performing numerical simulations of the 
system for ensembles ranging from N = 10^ to 10®. As customary with cou- 
pled phase oscillators, we rescale the constant k with the ensemble size as 
k — K/N, and fix X = 1. The frequencies Jlj of the internal phases are drawn 
from the distribution of Eq. (4). We fix fio = 0, and consider various values of 
the frequency dispersion a. Figure 8 shows snapshots of the ensemble taken at 
sufficiently long times, for three values of a. For small a (Fig. 8a), the phases 
(f)i are homogeneously distributed around the circle, and a high degree of corre- 
lation between 0j and the internal variables 6i is apparent. For large frequency 
dispersion, the distribution of phases is also approximately homogeneous but, 
on the other hand, the correlation with the internal variables is lost (Fig. 8c). 
At intermediate values of cr, we find strongly heterogeneous distributions for 
the phases (pi. Figure 8b shows a state where most of the ensemble is split into 
two clusters with opposite phases. 

To explain the appearance of these three regimes we first note that, for small 
0", the internal variables -and, therefore, the interactions- evolve slowly. Over 
such long time scales, the phases 0, are able to follow, almost adiabatically, 
the evolution of the interaction weights V{di, 9j). Numerical results show that, 
in this situation, the phase 0j and the internal variable 9i of each oscillator 
are related according to 

(Pi = ±0i + 00- (14) 
The sign factor of 9i and the value of 0o are the same for all the oscillators 
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Fig. 8. Long-time snapshots of an ensemble of 100 phase oscillators, for different 
values of the frequency dispersion: (a) a = 0.1, (b) a = 0.3, and (c) a = 0.5. 

in the ensemble. They are determined by the initial condition, and (po may 
slowly change with time. In the limits a ^ and N ^ oo, a homogeneous 

distribution of phases related to the internal variables as in Eq. (14) is a 
solution to Eq. (13). In fact, if such relation holds, we have 



For finite vahics of a, the relation between and 9i is maintained as long 
as the frequency dispersion remains small. For large a, on the other hand, 
most interaction weights show large changes over short times, and the phases 
are not able to adjust their values to such changes. In this situation, each 
oscillator is subject to rapidly fluctuating forces and, as a result, phases are 
homogeneously distributed, but no relation with the internal variables can 
persist. 

The transition between the ordered and the disordered regime, respectively 
found for small and large frequency dispersion, can be characterized by means 
of a suitably defined order parameter. With this aim, we first introduce for each 
oscillator the two-dimensional complex vector rrij = (exp[i(0i — 6*1)], exp[i(0i -|- 
Oi)]), and define the average 

m = ^^mi = (//+exp(iV'+),//-exp(iV'-)). (16) 



The average of the modulus // = |m| = + over sufficiently long times 
is an adequate order parameter. Indeed, if relation (14) is approximately ver- 
ified, we have /i. ~ 1. If, on the other hand, the values of 4>i and 9i are un- 
correlated, we find fx ~ 1/\/N. Figure 9 shows the numerical evaluation of 
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Fig. 9. The order parameter |U as a function of the frequency dispersion cj, for 
ensembles of various sizes N. The curve corresponds to the analytical evaluation of 
II, from Eq. (21). 

/i as a function of the frequency dispersion a, for different values of N. We 
note that, as N grows, the order parameter develops an abrupt inflection at 
(7 ~ 0.3, which suggests the presence of a critical phenomenon. As we show in 
the following, this critical phenomenon can be described analytically if a few 
assumptions, supported by numerical results, are made. 

Taking into account that, in the limit cr — > 0, the relation (14) is equiva- 
lent to the stationarity condition for an ensemble with homogeneous phase 
distribution, for larger values of a we can write 



where di{t) measures the deviation of each oscillator with respect to the sta- 
tionary state of small frequency dispersion. We assume that, for sufficiently 
large ensembles and long times, the distribution of these deviations is given 
by a well-defined function p{S), satisfying p{S) = p{—S). Implicitly, this also 
amounts to conjecture that the values of Si on the one hand, and 0j and 6i on 
the other, are not correlated. Replacing Eq. (17) into the definition of irij and 
using Eq. (16) to calculate the order parameter /x, we find 



Assuming that 0o in Eq. (17) is independent of time, Eqs. (13) and (18) yield 




(17) 




(18) 
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the following evolution equation for Sf. 

Si = fli — ^ sin Si- (19) 



This equation has been extensively discussed in connection with the dynamics 
of globally coupled phase oscillators [21]. For \fli\ < ///2, Si has two fixed 
points, at the solutions of 

sin Si ^ — -. 20 
1^ 



One of them is stable, and the other is unstable. On the other hand, for > 
n/2, there are no stable fixed points and \Si\ grows indefinitely with time. We 
can therefore distinguish between two subpopulations in the ensemble, with 
qualitatively different behaviour. Those oscillators for which \ili\ < n/2 reach, 
at asymptotically large times, a stationary deviation Si. In contrast, for those 
oscillator with > ///2 the phase 0, does not reach a stationary value with 
respect to the internal variable 9i. We refer to these two subpopulations as 
subensembles I and II, respectively. Whether a given oscillator belongs to one 
of the two subensembles depends on the value of According to Eq. (18), is 
given by the distribution of stationary deviations Si, which in turn depend on 
fjL through Eq. (20). Thus, the order parameter must be found self-consistently. 

The distribution of deviations p{S) consists of contributions from each sub- 
ensemble p — Pi+pii- Since the deviations Si for the oscillators of subensemble 
II are not stationary and move at different velocities, pii{S) is a flat distribution 
and, thus, does not contribute to fi in Eq. (18). The contribution of subensem- 
ble I, on the other hand, can be found from the identity pi{S)dS = g{Q)dQ, 
taking into account the connection between the deviation and the internal 
frequency of each oscillator, given through by Eq. (20). Replacing p{S) into 
Eq. (18) yields 

I g{n)Ji-^dn. (21) 

-n/2 ^ 



This is the self-consistency equation which makes it possible to flnd //., given 
the distribution of internal frequencies g{^). For our Gaussian distribution, 
Eq. (21) has a nontrivial solution > if cr < Uc = ^7r/32 « 0.313. As the 

frequency dispersion approaches the critical value 0"^ from below, the order pa- 
rameter vanishes as fi ~ |a — o-d^^^. For a > a^, the only solution is /i = 0. The 
curve in Fig. 9 shows the solution to Eq. (21) as a function of the frequency 
dispersion. The excellent agreement with numerical results validate our ana- 
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lytical approach. Remnant discrepancies can be ascribed to the enhancement 
of finite-size effects around the critical transition. 



Note that the same singular dependence of the order parameter close to 
the order-disorder transition would be found for a broad class of frequency 
distributions g{Q). In fact, since only small values of Q contribute to the 
integral in Eq. (21) when fi is small, the critical behaviour is determined by 
the shape of the maximum of g{^l) at ^2 = 0. Suppose that, for ?a 0, we 
have g{fl) ?a g[0) — 7^^. Replacing into Eq. (21) we find 




(22) 



Consequently, the critical point is determined by the equation g{0) = 4/tt. If, 
sufficiently close to the point where this equation holds, g{0) varies linearly 
with the control parameter a, the order parameter will behave as ~ \a — 
(7c I The critical exponent, here equal to 1/2, changes if the maximum of 
g{D,) at Q = is not quadratic. 

Using the complex vector m defined in Eq. (16), the equation of motion (13) 
for the phase of an individual oscillator can be recast as 

<i>i^^ sin(V'+ + - 0i) + ^ sin(V'- -e^-cPi). (23) 

The form of this equation emphasizes the mean-field nature of the interactions, 
as depends only on the individual variables (pi and 6i and on the averages 
involved in the definition of m. Since for asymptotically large ensembles, 
00, /X = + /i?. is zero for a > ac, both fx+ and fX- must also vanish beyond 
the transition. This implies that in the thermodynamical limit and for a > ac, 
— and the dynamics is thus frozen. For finite systems, //+ and /x_ are of 
order N'^^"^. The evolution of individual phases can therefore be thought as 
driven by random forces of that order. Below the transition, evolution becomes 
progressively slower as a approaches the critical point. 



3.3 Dynamical clustering 



As pointed out above, the dynamical regimes of small and large frequency 
dispersion are separated by a region, around a ~ 0.3, where clustering is 
found (Fig. 8b). To disclose the properties of this intermediate regime, we have 
numerically studied the evolution of the distribution of phases 0,. Recording 
the evolution of each phase, we construct a histogram where the height of 
the column with base (0, (p + A0) is proportional to the number of oscillators 
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Fig. 10. Density plot of the phase distribution (100-column histogram) as a function 
of time, for a system of iV = 10^ oscillators with frequency dispersion a = 0.3, along 
200 time units. Darker tones correspond to higher concentrations. 

whose phases he within that interval. Figure 10 shows a density plot of such 
histogram as a function of time, for a = 0.3. 

We note that a fraction of the ensemble is entrained into four clusters. A close 
inspection of these clusters shows that they are formed by the oscillators of 
subensemble II (see Sect. 3.2) with internal frequencies just above their lowest 
value, \fli\ > ii/2. Though the clusters are not sharply localized, it is still 
possible to define their phase by averaging (pi over the oscillators belonging 
to each cluster. Two of these clusters, whose phases differ by tt, move at a 
well defined velocity in one direction. The other two, whose phases also differ 
by TT, move in the opposite direction and with the same speed. The two pairs 
of clusters are formed by oscillators whose internal frequencies have opposite 
signs. Figure 11a illustrates the situation. Due to their contrary motion, the 
two pairs of clusters cross each other recurrently. When these crossings happen, 
their populations are temporarily superimposed in phase. Consequently, two 
big clusters with opposite phases are recurrently built up (Figs. 10 and lib). 
Figure 8b shows a 100-oscillator ensemble at one of such events. 

The recursive formation of the two anti-phase big clusters can be quanti- 
tatively characterized by an order parameter. The modulus of the complex 
quantity 

^W = 4Eexp(2i0,) (24) 
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Fig. 11. (a) Schematic representation of the relative position and motion of the two 
pair of clusters found in the clustering regime, (b) Histogram of the distribution 
of phases when the two cluster pairs cross each other, in a system of N = 10^ 
oscillators with frequency dispersion a = 0.3. 

is approximately equal to one if the whole ensemble splits into two anti-phase 
well-localized clusters of similar sizes, while for incoherent or higher-order 
clustered states we have \z\ ~ N"^/"^. Figure 12a shows the evolution of 
in the clustering regime, displaying its irregular oscillations between small 
and large values. For the same realization. Fig. 12b shows the phase difference 
= (pa — 4>b for two oscillators a and b in subensemblc II. Large positive 
values of cos A0 identify states where the two oscillators are part of the same 
cluster, while negative values correspond to the situation where they are in 
different clusters. The intermittent transitions of cos Acj) between positive and 
negative values demonstrate the dynamical nature of clustering, where clusters 
maintain their identity but elements can migrate between them. 

The average of \z{t)\ over long times provides a suitable order parameter, C, 
for the detection of clustered states. Note that ( and are independent order 
parameters, as they measure Fourier contributions of different order to the 
distribution of phases. Figure 13 shows numerical results for ( as a function 
of the frequency dispersion. As a approaches the critical value cXc from below, 
C grows and exhibits a peak just at the left of the critical point. The profile 
of this peak becomes better defined as the ensemble size grows. The order 
parameter grows from zero at o" ~ 0.2 and reaches the maximum value ( ^ 0.25 
for a < (Je. Beyond the transition, ( drops abruptly. Results for the largest 
ensembles suggest that, in the thermodynamical limit, it vanishes for a > Of.- 
The suppression of clustering beyond the transition is a consequence of the fact 
that, as discussed in Sect. 3.2, the evolution of asymptotically large ensembles 
becomes frozen. As a approaches the critical point and the order parameter [i 
tends to vanish, the typical evolution time scales -which, according to Eq. (23), 
are of order are increasingly large. Thus, the formation of clusters is 

progressively retarded and, finally, it is suppressed. 
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Fig. 12. (a) Evolution of \z{t)\, Eq. (24), for a system of = 10^ oscillators wit h 
frequency dispersion a = 0.3. (b) Phase difference of two oscillators of subensemble 
II with similar frequencies, in the same realization as in (a). 
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Fig. 13. The order parameter as a function of the frequency dispersion a for 
oscillator ensembles of various sizes. The vertical dashed line indicates the critical 
point of the order-disorder transition, ctc = y^Tr/S^. 

The appearance of the four clusters in the transition regime can be quali- 
tatively understood taking into account that, for an oscillator with internal 

frequency Qi, the most persistent interactions affecting its dynamics are due 
to oscillators with similar frequencies, Qj ^ fij. In such case, in fact, the in- 
teraction weight V{9i, 9j) remains almost constant over very long time scales, 
of order — Fixing the attention on a group of oscillators with similar 
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frequencies, we realize that their interactions arc weighted by almost time- 
independent factors which, however, can take a different value for each oscil- 
lator pair. We can have both positive and negative interaction weights. Now, 
it is known that an ensemble of coupled oscillators with disordered interac- 
tions forms, if the level of frustration is moderate, two anti-phase localized 
clusters [18]. In our system, triggering the formation of these two clusters re- 
quires moreover that the phases </)j are not pinned to the internal variables 
9i, and that a sufficiently large population with similar frequencies actually 
exist. These two conditions are met by the oscillators of subensemble II with 
frequencies close to the hmit \Q\ — /i/2, where the density is larger. We ex- 
pect that two clusters form in the region of positive frequencies, and other two 
clusters appear for i7 < 0. The symmetry of this configuration insures that the 
two pairs of clusters will move in opposite direction and with the same speed. 
As far as a remains below the transition point, increasing the frequency dis- 
persion contributes in two cooperating ways to the growth of the population 
of these clusters: larger values of a imply more oscillators with higher internal 
frequencies and, at the same time, a larger fraction in subensemble II, as /i 
decreases. 



4 Conclusion 

We have studied, both numerically and analytically, the dynamical regimes of 
an ensemble of coupled motile elements whose interactions are determined by 
the evolution of their internal states. The interaction of each pair of elements is 
alternatively attractive or repulsive, changing its sign within time scales asso- 
ciated with the relative variation of the respective internal variables. Though 
in our study we have considered specific forms both for the spatial dependence 
of the interactions and for their modulation by the internal state, our main 
results are expected to hold for a broad class of similar models. The following 
summary is focused on such generic conclusions. 

Transitions between different dynamical regimes are controlled by two param- 
eters: the range of interactions and the typical evolution time of the internal 
variables. The latter determines, in turn, the time scales within which inter- 
actions change their sign. The simultaneous presence of both positive and 
negative pair interactions induces the formation of localized clusters of mu- 
tually attracting elements. The remnant repulsion drives clusters away from 
each other, a tendency to expansion that must be counteracted by suitable 
boundary conditions. Due to the change of the sign of interactions, cluster- 
ing is dynamical, as elements intermittently aggregate, disperse, and migrate 
between clusters. The actual development of localized clusters, however, re- 
quires that the interaction range is neither too small nor too large. Clusters 
are practically not formed if that range is smaller than the average distance 
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between elements, while for large interaction ranges other kinds of spatial or- 
ganization, such as "strings" of elements and condensation in the boundaries, 
become more frequent. 



When the range of interactions is large as compared with the spatial size of the 
system, interactions act in conditions similar to those of global coupling [17], 
though different interaction weights are still possible between different element 
pairs. In such conditions, the spatial distribution of the ensemble becomes less 
relevant to the collective dynamics, because interactions are now independent 
of the relative position of elements. Yet, their individual positions are still 
determined by their interactions, which are now mainly driven by the internal 
variables. When the evolution of the internal variables is sufficiently slow, the 
spatial distribution can follow it almost adiabatically, and the relative position 
of the elements adjusts to their instantaneous internal state. As a result, a 
state of high correlation between spatial and internal variables develops. Due 
to the autonomous internal evolution of each element, this ordered state is 
dynamical, and the spatial distribution keeps varying with time. On the other 
hand, if the internal variables change too fast, the spatial coordinates fail 
to follow their evolution, and correlation between them breaks down. If the 
internal evolution is sufficiently fast and the time average of pair interactions 
is zero, the spatial distribution of the ensemble becomes frozen. In our specific 
model, the transition between the ordered and the disordered state has the 
properties of a critical phenomenon. In the disordered state, evohition is frozen 
in the limit of an infinitely large ensemble. The transition is mediated by a 
regime of dynamical clustering, induced by the persistent interaction between 
elements whose internal variables change within similar time scales, and whose 
positions and internal states are not strongly correlated. 



The dynamics of internal states is expected to be an essential ingredient in 
the description of systems whose collective evolution resembles the complex 
collective dynamics of biological and social populations [14]. Spontaneous de- 
velopment of coherent behaviour may be enhanced if individual internal vari- 
ables are allowed to interact with each other. In fact, the organization of the 
ensemble at the level of the individual internal states, due to their mutual 
interaction, would imply coherence also at the level of other variables, such as 
the spatial position. These processes are expected to be relevant to the mod- 
eling of social-like mechanisms, of the type of imitation, social influence, and 
formation of coalitions and hierarchical structures. The investigation of their 
role in the emergence of order has already been initiated in the literature [15], 
also in connection with evolutionary adaptive phenomena [22,23,24]. 
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